Immune Infiltration Correlates with Transcriptomic Subtypes in Primary ER+ Invasive Lobular Breast Cancer

Abstract Understanding interplay of breast cancer and microenvironment is critical. Here, we identified two transcriptomic subtypes and five immune infiltration patterns from RNA-seq and multiplex immunohistochemistry from 21 ER+/HER2- invasive lobular breast cancers. The proliferative subtype associated with increased immune infiltration especially by immunosuppressive regulatory T-cells and macrophages. We also defined a TAM-Low signature, which associated with lower infiltration of proliferative, pro-inflammatory TAM, and improved outcome in patients with ER+ tumors.


Main text
Immunotherapy (IO) is becoming an important option in breast cancer treatment.In triplenegative breast cancer (TNBC), pembrolizumab showed potential benefits in high-risk earlystage disease and metastases (1,2).In estrogen receptor-positive, HER2-negative (ER+/HER2-) early-stage breast cancer, tumors with lower ER and higher PD-1 were more likely to respond to neoadjuvant immune checkpoint inhibitors (ICI) (3,4).However, most trials excluded invasive lobular breast cancer, and the investigation of clinical biomarkers beyond ER, PD-1, or TILs is limited.Transcriptomic data is useful to identify prognostic biomarkers, including tumor subtypes and immune patterns.The intrinsic molecular subtype was one of the first classification systems derived from patient tumor RNA expression (5).In TCGA, three ILC subtypes were revealed from tumor RNA-seq data by unsupervised clustering, where the reactive-like subtype had better disease free survival and overall survival than the proliferative subtype (6).Similar classifications were also derived in the RATHER study (7).More recently, four transcriptomic subtypes were defined in HR+/HER2-breast cancer, among which the receptor tyrosine kinase (RTK)-driven subtype had poorest outcomes compared to the canonical luminal, immunogenic, and proliferation subtypes (8).For immune patterns, immune quantification data is scarce, therefore, to relate to clinical prognosis, studies often calculated scores from marker gene expression of individual immune cell types in bulk RNA-seq data.For example, one study reported that higher immune infiltration was linked with lower tumor purity, and predicted better overall survival (9).Another study analyzed TCGA and METABRIC data, and found that tumors with higher M1-like signature or M1/M2 ratio were likely more aggressive with poorer survival (10).However, prognostic impact of immune infiltration patterns remains ambiguous, considering decreased estimation accuracy for more complex immune cell subtype composition.Furthermore, most studies excluded ILCs, which had distinct immune infiltration than no special type (NST) tumors (11,12), thus introducing bias especially in data-driven conclusions.Further research is thus warranted to identify prognostic immune infiltration patterns in ILC tumors specifically.In this study, we performed integrative analysis of both RNA-seq and multispectral immunohistochemistry (mIHC) data in primary ER+/HER2-ILC.We identified two transcriptomic ILC subtypes, five immune infiltration patterns, and linked the ILC subtypes with immune infiltration.We then correlated gene expression with immune cell density, and defined a tumor associated macrophage (TAM)-Low signature, which was negatively associated with TAM infiltration.We further showed in public datasets that TAM-Low was inversely correlated with presence of proliferating TAMs, and predicted improved survival outcomes in ER+ breast cancer.21 primary ER+ treatment-naïve lobular breast tumors were included in our cohort (Supplementary Table 1).We performed RNA-seq in all samples, and mIHC in 13 of the tumors using the same FFPE block (Fig. 1A).We identified CDH1 truncating mutations in 10 samples from RNA-seq (Supplementary Table 2, Supplementary Figure 1A).We then computed transcriptomic subtypes.Using TCGA classifier (6), we identified 12 proliferative, three reactivelike, and six immune-related tumors, and combined the latter two groups as the non-proliferative subtype (Fig. 1B).Consistently, we found higher expression of cell-cycle pathways in the proliferative group, and higher expression of immune-related pathways in the non-proliferative subtype, using reference pathways from the RATHER ILC cohort (enriched pathways of hormone-related and immune-related subtypes, respectively) (Fig. 1C).Specifically, the nonproliferative subtype had higher stromal score and mostly consisted of luminal A and normal-like subtypes.In contrast, the proliferative subtype had higher tumor cell purity, lower stromal composition, and was composed of multiple PAM50 types, including luminal A, luminal B, basal, and HER2-enriched subtypes (Fig. 1D-E).Next, we examined mIHC immune infiltration patterns from ROIs in tumor and stromal regions (Fig. 2A).In general, macrophage showed highest infiltration, followed by CD4 T cells and CD8 T cells, all with higher infiltration in stromal than tumor regions (Fig. 2B).Correlation based on infiltration density of the different immune cells, reflecting their associations within or between tumor and stroma, we observed three clusters -stromal and tumor B cells, Treg cells and macrophages, and CD4 T cells and CD8 T cells (Fig. 2C).Furthermore, we identified five immune infiltration patterns in both stromal and tumor regions from non-negative matrix factorization (NMF) clustering (Fig. 2D).In stromal regions, pattern 1, 3, and 5 had generally low infiltration of all immune cells, pattern 2 had higher macrophages and Tregs, and pattern 4 had higher infiltration of CD4 T cells and CD8 T cells (Fig. 2E-F).Similar immune infiltrating patterns were observed in tumor regions (Supplementary Figure 1D-E).This suggested a potential immunosuppressive microenvironment in pattern 2 while more immune active environment in patterns 4. Specifically, within each patient, the tumors had up to 5 immune patterns, and was often dominated by 1 (blue arrows) or 2 (orange arrows) co-existing patterns (Supplementary Figure 1B).We then linked immune infiltration with the ILC transcriptomic subtypes, defined in the previous section.The proliferative subtype showed higher infiltration of almost all immune cell types at both stromal and tumor areas than the non-proliferative subtype (Fig. 2G-H).Specifically, pattern 2 (likely immunosuppressive) predominated in the proliferative subtype (Supplementary Figure 1C).Altogether, it suggested higher immune infiltration, and potentially more immunosuppressive cell types in the proliferative than non-proliferative ILCs.Access to both mIHC as well as mRNA expression data allowed us to analyze correlations between gene expression and immune cell infiltration.We therefore calculated Spearman's correlation between expression of each individual gene (26,485 genes, TPM) and infiltration of major immune cell type (5 types) in stromal or tumor regions (median number per mm 2 among ROIs), including all possible pairwise combinations.Filtering by FDR<0.05,we identified 651 significant correlations between expression of genes and immune cell subtypes (Fig. 3A).The majority of correlations were negative associations between gene expression and macrophage infiltration.Combining genes inversely related with presence of stromal or tumor macrophages, we generated a 'TAM-Low' signature of 483 genes (Fig. 3A, Supplementary Figure 2A), which was enriched in extracellular matrix pathways (Fig. 3B).Next we compared signature expression of TAM-Low with marker genes of multiple subtypes of tumor-associated macrophages (Supplementary Table 3) from literature in ER+ primary breast cancer, using RNA-seq of METABRIC and SCAN-B (13)(14)(15).In both datasets, the TAM-Low signature was negatively correlated with proliferating TAMs (Prolif-TAMs) (Fig 3C ), a proinflammatory subtype characterized by high expression of HMGB1 and cell cycle-related genes (CCNA2, CDC45, CDK1, MKI67, RRM2, STMN1, TOP2A, and TYMS).In summary, we identified a negative association between the TAM-low signature and tumor associated macrophages, and specifically our data suggested that TAM-low genes such as those from the extracellular matrix (ECM) pathway might hinder the infiltration of proliferating, pro-inflammatory tumor associated macrophages.And finally, we investigated association of TAM-Low with both survival and therapy response.In ER+ breast cancer, the TAM-Low signature predicted improved relapse free survival as well as overall survival in METABRIC ER+ tumors, adjusting for multiple demographic and clinical variables in multivariate Cox regression (Fig. 3D-E).Result in SCAN-B ER+ tumors was similar but lacked significance (Supplementary Figure 2C).For therapy response, we analyzed the TAM-Low signature levels in ER+/HER2-primary breast tumors from patients treated with neoadjuvant aromatase inhibitor (AI) in the POETIC trial (16).Tumors from AI responders had higher median value TAM-Low signature expression than non-responders (Fig. 3F-H) at baseline and more pronounced after neoadjuvant retreatment, although this did not reach significance.In line with our prior observation of a negative association between the TAM-Low signature and marker genes of proliferating TAMs (Prolif-TAMs), we observed significantly lower proliferating TAMs (Prolif-TAM) after AI treatment in the responders.Using expression of marker genes of other TAM subtypes (15), we failed to identify another set of TAMs significantly different between non-responders and responders, except for lower interferon-primed TAMs (IFN-TAM) at baseline in responders (Fig. 3F).In summary, this data suggests that the TAM-low signature is associated with improved endocrine therapy response and outcome ER+ breast cancer.Finding biomarkers is essential to guide precision immune therapy in breast cancer.In TNBC, higher tumor-infiltrating lymphocytes and mutation burden predicted higher ICI response rate.In early HR+/HER2-cancer, which had relatively lower mutation incidence and immune infiltration, positive lymph node, ER%<10%, and combined positive score 1 were associated with higher pCR rate with neoadjuvant pembrolizumab treatment (3), and elevated PD-L1 scores, ER%<50%, PR%<10%, and stromal TILs 1% indicated pCR and no residual cancer burden with neoadjuvant nivolumab therapy (4).However, clinical associations of other immune cell types or genes were understudied, and whether pCR translates to long-term survival benefit remained unclear.Moreover, clinical trials did not focus on lobular breast cancer, which had a somewhat distinct immune landscape than NSTs (11,17,18).Our study focused on primary ER+/HER2-early-stage ILC, and found ILC transcriptomic subtypes were closely related with immune infiltration phenotypes.Specifically, in the proliferative ILC subtype greater infiltration of multiple immune cell types, higher occurrence of immuno-suppressive patterns was also observed.This was consistent with clinical observations that high TIL level in ILC was associated with young age, proliferative tumors, and lymph node invasion (18).From our previous study, higher macrophage infiltration was related to lower recurrence rate in IDC, but this was not observed in ILC (11).In general, large-scale data is lacking for clinical outcome associations with TIL subtypes or other immune cell types.This hence warrants future effort to incorporate mIHC quantification in prospective cohorts, and to include ILC cases in clinical trials (17).Most current immunotherapies focus on lymphocytes.However, macrophages, rather than T cells, were dominant in ILCs (11).In our study, we also found macrophage infiltration as the 'hub' in association with tumor gene expression.Our TAM-Low signature which was directly derived from correlation analysis between gene expression and immune infiltration based on mIHC has potential as a promising biomarker, predicting low infiltration of proinflammatory Prolif-TAMs, and improved survival in patients with ER+ breast cancer.Furthermore, higher Prolif-TAMs likely marks intrinsic aromatase inhibitor resistance in ER+/HER2-primary breast tumors.As resident tumor macrophages expressing MKI67, CDK1, and CDC45, HMGB1, and others, Prolif-TAMs were likely pro-inflammatory, profibrotic, and promoting tumor progression (15).In pancreatic ductal adenocarcinoma, targeting Prolif-TAM promoted anti-tumor response via cytotoxic CD8+ T cell redistribution (19).Interestingly, from our original IDC and ILC cohort, we had reported a subgroup of cycling monocytes/macrophages (cluster 4) expressing MKI67 and SPP1 genes, which likely corresponded to Prolif-TAMs (11).Of note, each TAM subtype affects breast cancer differently.For example, the FOLR2+ tumor-stroma-resident macrophages correlated to higher T cell infiltration and better prognosis (20).Conversely, TREM2+ macrophages likely determined an immunosuppressive metastatic niche, with TREM2 predicting poorer overall survival and relapse-free survival in triple-negative breast cancer (21,22).Therefore, we encourage quantification of individual TAM subtypes in breast cancer and specifically in ILC as a more accurate immune phenotyping in clinical response prediction.Limitations of this study include small cohort size, variations in ROI selection, and absence of immunophenotyping panel of a broader range of immune subtypes, including the macrophage population.There are also limitations to combining different methodologies for analysis of immune cell infiltration and activities, as reflected by some discrepancy between mIHC and RNA-seq deconvolution algorithms -e.g., proliferative ILCs had higher stromal or immune scores and showed lower expression of immune signaling than non-proliferative ones, including the JAK-STAT pathway and local acute inflammatory response (LAIR pathway).A possible explanation is that immune signals were expressed by tumor cells, thus leading to spillover effect and elevated immune scores in non-proliferative tumor (23).
In conclusion, through integrated RNA-seq and mIHC analysis of ER+/HER2-primary ILC tumors, we found that the proliferative ILC transcriptomic subtype exhibited higher immune infiltration.We would like to propose the TAM-Low signature as a potential novel biomarker for predicting reduced infiltration of pro-inflammatory and pro-tumorigenic macrophages and improved survival in ER+ breast cancer, and encourage its further analysis in additional studies.Future studies are warranted to quantify TAM subtypes in ILC, as indicators of endocrine treatment response and survival, as well as to improve potential use of immunotherapy.

Study design, sample acquisition, and sample processing
Patient recruitment was described previously (11).For the 21 cases, formalin-fixed paraffinembedded (FFPE) blocks of primary ER+/HER2-ILC tumors were acquired from University of Pittsburgh Biospecimen Core, and clinical information was queried from the UPMC Cancer Registry.For each block, three sections (at the top, bottom and middle of the block) were stained with hematoxylin and eosin (H&E) and inspected by pathologist to circle out regions with tumor epithelial cell proportion 40%.5 sections (of 10uM thickness each) were used for RNA extraction and sequencing at these circled regions.And another 3 sections (of 4uM thickness each) were used for multispectral immunophenotyping by mIHC.Protocol of this study was approved by Institutional Review Boards at University of Pittsburgh (protocol numbers PRO17100602 and PRO16060685).

RNA sequencing and bioinformatics analysis
Sample RNA was extracted with Qiagen tissue DNA/RNA kit from FFPE sections, assessed with High Sensitivity RNA TapeStation for quality and concentration, and sequenced with Illumina Nextseq 500 75 bp pair-end in 150 cycles at University of Pittsburgh Health Sciences Sequencing Core.For RNA-seq data, we used FastQC for FASTQ data quality control, trimmomatic for adapter trimming, Hisat2 for alignment to human genome 38 (hg38), and HTSeq for quantification, to generate normalized expression in transcripts per million (TPM).To call CDH1 mutations from RNA-seq data, we used GATK v4 MuTect2 pipeline (broadinstitute.org) to call all somatic SNVs and indels from bam files calculated with STAR (v2.7.5a), and selected CDH1 mutations from annotation from ANNOVAR (24).For ILC subtyping, we first trained a Nearest Centroid classifier in original TCGA ILC samples with RNA expression of the 60-gene signature, using original annotation of the proliferative, reactive-like, and immune-related subtypes (6).We then applied the trained classifier to RNA-seq data of this study, which generated one subtype for each sample.Both training and testing sets were normalized in log2TPM and transformed into gene-wise z-score.Finally, we combined the reactive-like and immune-related subgroups as the 'non-proliferative' subtype, due to limited sample numbers of both classes.We used GSVA (v1.48.2, ssGSEA) to calculate signature scores of immune-related and cell-cycle related pathways from TPM of our samples.We calculated PAM50 subtype probability with Genefu (v2.26.0), and estimated tumor purity, stromal score, and immune scores using deconvolute_estimate (immunedeconv v2.1.0).

mIHC staining, quantification, and statistical analysis
Workflow of mIHC was described with processed data (median infiltration of each immune cell type at patient level) previously (11).In brief, FFPE sections of 4uM thickness were stained with an immunophenotyping panel (CD20, CD4, CD8, FOXP3, CD68, Pan-CK), from which immune cells were counted using InForm software from randomly selected regions of interests (ROI) (3-13 / section) at representative immune infiltration areas, from Vectra Polaris (PerkinElmer) scans at x10.Specifically, infiltration of each immune cell type was calculated in tumor and stromal regions respectively (number/mm 2 ).We computed Spearman's correlation for abundance of each pair of immune cell types (number/mm 2 ) and used FDR<0.05 to select statistically significant correlations.Non-negative matrix factorization was used to cluster ROIs based on Spearman's correlation matrix of log2 normalized counts (BigNmf, v1.0.0).from which 5 clusters were selected based on elbow point of cophenetic correlation coefficient.The corresponding consensus matrix was plotted in Euclidean distance in Ward's linkage.

TAM-Low signature generation and pathway analysis
We calculated Spearman's correlation between expression of each gene (TPM) and abundance of each immune cell type (number/mm 2 ) at tumor and stromal region (median values among ROIs).A gene-immune cell pair correlation is significant if with FDR<0.05.Genes which were significantly negatively correlated with stromal or tumor macrophage infiltration were included into 'TAM-Low' signature (483 genes, Fig. 3, Supplementary Figure 2).Pathway analysis was performed with EnrichR, using Gene Ontology, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome datasets (25).

Figure 1 .
Figure 1.Study design and ILC subtype identification from RNA-seq A. Study design and workflow.Treatment-naïve ER+/HER2-primary ILC tumors were processed as FFPE blocks, from which FFPE sections underwent RNA-seq in tumor-rich regions selected by pathologist (N=21), and multispectral IHC (mIHC) with immunophenotyping panel (N=13).For mIHC, immune cell infiltration was quantified in regions of interest (ROIs) in both tumor and stromal area.B. Expression of ILC classifier genes in the ER+/HER2-ILC cohort, showing the two ILC subtypes (proliferative in orange box, non-proliferative (immune + reactive) in blue box), PAM50 subtypes, and other demographic and clinical information.C. Signature score (ssGSEA score) of reference ILC subtype pathways from RATHER study.Only pathways with significant signature score difference between non-proliferative and proliferative ILC subtypes were shown (Mann-Whitney-Wilcoxon test).D. Estimated tumor purity, stromal score, and immune score from PUREE and ESTIMATE algorithm.E. PAM50 subtype composition of non-proliferative and proliferative ILC tumors.

Figure 2 .
Figure 2. Immune cell infiltration landscape from mIHC A. Immune cell infiltration in tumor and stromal regions among ROIs (number/mm 2 ), annotated with ILC cluster, PAM50 subtypes, and demographic or clinical information.B. Histogram of immune cell type abundance (/mm 2 ) in stromal and tumor region ROIs.C. Spearman correlation (FDR<0.05) of immune cell infiltration (number/mm 2 ) between cell types.D. NMF consensus matrix from immune cell infiltration in both stromal and tumor regions among ROIs, showing 5 immune infiltration patterns.E. Stromal region immune cell infiltration (number/mm 2 ) among immune infiltration patterns, showing each pattern per panel.F. Stromal region immune cell infiltration (number/mm 2 ) among immune infiltration patterns, showing each immune cell type per panel.****p<0.0001,***p<0.001,**p<0.01,*p<0.05,comparison (horizontal brackets) with no label, p0.05.G-H.Immune cell infiltration (number/mm 2 ) in proliferative and non-proliferative ILC tumors in stromal (G) and tumor regions (H).

Figure 3 .
Figure 3. Derivation and clinical associations of TAM-Low signature A. Number of gene-cell pairs with significant (FDR<0.05)positive or negative Spearman's correlation between gene expression (TPM) and immune cell infiltration (median value among ROIs in number/mm2).B. Pathway enrichment of TAM-Low signature genes (483 genes).C. Spearman correlation (FDR<0.05) between TAM-Low and other TAM subtype signatures in Figure 1